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Abstract 

A  key  ingredient  in  simulation  of  flow  in  porous  media  is  accurate  determination  of  the 
velocities  that  drive  the  flow.  Large-scale  irregularities  of  the  geology  (faults,  fractures,  and 
layers)  suggest  the  use  of  irregular  grids  in  simulation.  This  paper  presents  a  control- volume 
mixed  finite  element  method  that  provides  a  simple,  systematic,  easily  implemented  proce¬ 
dure  for  obtaining  accurate  velocity  approximations  on  irregular  block-centered  grids.  The 
control- volume  formulation  of  Darcy’s  law  can  be  viewed  as  a  discretization  into  element¬ 
sized  ''tanks”  with  imposed  pressures  at  the  ends,  giving  a  local  discrete  Darcy  law  analogous 
to  the  block-by-block  conservation  in  the  usual  mixed  discretization  of  the  mass-conservation 
equation.  Numerical  results  in  two  dimensions  show  second-order  convergence  in  the  veloc¬ 
ity,  even  with  discontinuous  anisotropic  permeability  on  an  irregular  grid.  The  method 
extends  readily  to  three  dimensions. 
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Contract  No.  NAS  1-1 9480  while  the  first  author  was  in  residence  at  the  Institute  for  Computer  Applications 
in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001 
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Figure  1:  Logically  rectangular  grid  of  irregular  quadrilaterals 


1  Introduction 

As  techniques  of  reservoir  description  become  more  sophisticated,  it  becomes  increasingly 
important  to  model  flows  of  reservoir  fluids  accurately.  In  particular,  it  is  desirable  to 
accurately  represent  large-scale  irregularities  of  reservoir  geology  in  models.  The  control- 
volume  mixed  finite  element  method  allows  the  use  of  irregular  grids  while  maintaining 
many  of  the  familiar  properties  of  block-centered  finite  difference  methods  for  rectangular 
grids.  For  example,  it  preserves  the  notion  of  block-by-block  material  balance,  with  physical 
interblock-flow  terms.  It  also  yields  an  analogue  of  the  local  discrete  Darcy  law  (relating  a 
combination  of  fluxes  to  a  pressure  drop  between  blocks)  on  a  block-sized  ‘‘tank”  between 
two  pressure  nodes.  The  method  can  be  applied  to  any  “logically  rectangular”  grid  of 
irregular  quadrilaterals  in  two  dimensions,  or  analogous  hexalaterals  in  three  dimensions. 
In  two  dimensions,  “logically  rectangular”  means  that  each  grid  block  can  be  assigned  an 
index  {i,j)  such  that  it  shares  an  edge  with  the  usual  and  {i,j  ±  1).  An  example  of 

such  a  grid  is  shown  in  Figure  1.  Even  on  rectangular  grids,  the  control- volume  mixed  finite 
element  method  can  be  more  accurate  than  finite  differences.  However,  its  main  advantage 
is  that  it  can  be  applied  in  a  simple,  straightforward  way  to  obtain  accurate  results  on 
irregular  grids,  allowing  the  permeability  coefficient  to  be  anisotropic  and  discontinuous. 
Formulations  in  common  use  in  the  petroleum  industry,  such  as  corner-point  geometry  [18], 
are  well-known  to  be  inconsistent  on  general  logically  rectangular  grids.  Since  the  motivation 
for  irregular  grids  is  the  accuracy  of  the  discretization,  such  formulations  are  best  limited 
to  careful,  restricted  use  in  practice. 

We  describe  the  method  for  rectangular  grids  in  Section  2,  and  for  irregular  quadrilateral 
grids  in  Section  3.  The  method  is  presented  in  the  context  of  a  model  pressure  equation. 
Section  4  contains  results  from  numerical  experiments,  including  a  comparison  to  block- 
centered  finite  differences  and  numerical  convergence  results.  Section  5  is  a  summary. 

We  conclude  this  introductory  section  with  brief  descriptions  of  mixed  and  control- 
volume  finite  element  methods,  two  methods  that  provide  some  of  the  building  blocks  for  the 


control- volume  mixed  finite  element  method.  Between  these  descriptions  we  also  summarize 
some  other  methods  recently  proposed  for  irregular  grids  in  porous  media.  These  other 
methods  share  our  goal  of  circumventing  the  inconsistencies  of  approaches  such  as  corner- 
point  geometry. 

Mixed  finite  element  methods.  We  recall  some  of  the  essentials  of  mixed  finite- 
element  methods.  The  idea  is  to  represent  a  partial  differential  equation  as  a  system  of 
lower-order  equations,  solving  these  for  multiple  variables  of  physical  interest.  To  keep  the 
description  simple,  assume  incompressible  single-phase  flow,  neglecting  gravitational  effects, 
so  that  the  pressure  equation  takes  the  form 

-V-(^Vp)  =  g,  x€fi,  (1) 

where  k  (scalar  or  anisotropic  tensor)  is  the  permeability,  /i  the  fluid  viscosity,  p  the  pressure, 
q  a  source/sink  (e.g.,  well)  term,  and  Q  is  the  reservoir  with  boundary  dO,.  Also  for  simplicity, 
take  the  no-flow  boundary  condition 


Vp-n-0,  xedQ.  (2) 

A  mixed  method  separates  Darcy’s  law, 

V  =  -tvp,  (.3) 

where  v  is  the  velocity  vector,  from  conservation  of  ma.ss, 


V-v  =  9,  (4) 

and  solves  the  system  Eqs.  (3)-(4)  for  v  and  p,  instead  of  solving  Eq.  (1)  for  p  and  applying 
Eq.  (3)  to  obtain  v. 

For  the  standard  (not  control-volume)  mixed  method,  following  [19],  write  Eq.  (3)  in  the 
form  {p/k)v  +  'Vp  =  0,  multiply  by  a  vector  test  function  w,  integrate  over  fi,  and  integrate 
by  parts  to  obtain 

J  •  wdx  -  y  V  •  wpdx  =  0.  (5) 

Multiply  Eq.  (4)  by  a  scalar  test  function  2  and  integrate  over  to  see  that 


/  V  V2dx  =  f  qz 
Jq  Jn 


dx. 


(6) 


The  requirements  of  w  and  z  are  that  z  and  the  components  of  w  be  square-integrable,  that 
the  divergence  V  •  w  be  square-integrable,  and  that  w  •  n  =  0  on  the  boundary  SQ.  At  this 
point,  the  dififerentiai  equations  are  still  being  viewed  continuously,  with  p  and  v  satisfying 
the  same  conditions  as  z  and  w,  respectively. 

In  two  dimensions,  the  discrete  Raviart-Thomas  elements  can  be  rectangles  or  triangles. 
In  either  case,  for  the  lowest-order  elements  that  are  analogous  to  block- centered  finite 
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Figure  2:  Velocity  basis  functions  on  rectangles  and  triangles 


differences,  p  and  2:  are  piecewise-constant.  The  velocity  functions  v  and  w  can  best  be 
viewed  by  associating  a  degree  of  freedom  with  the  flux  (normal  component  times  edge 
length)  on  each  inter-block  edge;  this  covers  both  the  rectangles  and  the  triangles.  Velocity 
functions  for  which  the  flux  is  1  on  one  edge  (hence  the  normal  component  is  l/l^'l,  where 
\E\  is  the  length  of  the  edge)  and  0  on  other  edges  are  pictured  in  Figure  2.  Fluxes  vary 
linearly  in  the  direction  of  the  velocity.  Continuity  of  flux  is  assured  in  either  case. 

There  are  two  principal  advantages  to  this  approach.  First,  with  piecewise-constant  z, 
Eq.  (6)  yields  conservation  of  mass  on  an  element- by-element  basis,  in  analogy  with  block- 
centered  finite  differences.  Second,  approximating  v  directly  by  finite  elements  can  be  much 
more  accurate  than  solving  for  p  and  invoking  Eq.  (3),  especially  when  the  mobility  k/ p  is 
not  smooth.  An  example  of  the  importance  of  this  appears  in  [11],  where  a  mixed  method 
avoided  spurious  viscous  fingering  effects  that  had  appeared  in  Galer kin  finite-element  results 
for  miscible  displacement.  On  standard  grids,  a  mixed  method  has  also  reduced  numerical 
dispersion  in  an  industry  simulator  [10].  For  standard  grids,  the  convergence  and  accuracy  of 
these  mixed  methods  are  well-established,  both  independently  [19]  and  as  part  of  a  coupled 
system  for  miscible  displacement  [15].  This  theory  extends  in  a  straightforward  way  to  a 
grid  of  parallelograms,  which  are  linear  images  of  rectangles. 

For  arbitrary  quadrilateral  grids,  Raviart  and  Thomas  [19]  and  Thomas  [23]  defined 
appropriate  pressure  and  velocity  spaces  via  the  Piola  mapping  (see  Section  3.1).  Stability 
and  convergence  were  demonstrated  more  recently  by  Wang  and  Mathew  [25]  and  Shen  [22]. 
These  error  estimates  are  of  significant  interest  mathematically,  but  of  limited  practical 
value  in  the  geological  context  motivating  our  work,  because  they  depend  on  continuity 
of  the  normal  and  tangential  components  of  v  at  an  interface;  at  interfaces  where  k  is 
discontinuous,  the  tangential  component  of  v  will  also  be  discontinuous.  Farmer  ei  ai  [12] 
have  applied  these  methods  to  petroleum  reservoir  simulation. 

Other  methods  for  irregular  grids.  In  the  petroleum  industry,  Aavatsmark  et  al.  [1, 
2]  have  proposed  and  tested  a  block-centered  finite-difference  method  involving  partitioned 
fluxes  between  subdivisions  of  irregular  blocks.  The  objective  is  continuity  of  pressure 
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and  velocity  with  suitable  interpolation  between  nodes;  this  would  be  overdetermined,  so 
some  constraints  are  relaxed.  Of  a  similar  flavor  is  a  method  proposed  and  analyzed  by 
Thomas  and  Trujillo  [24],  which  is  a  mixed  discretization  with  Eq.  (3)  written  in  the  form 
V  4-  (A:///)Vp  ==  0  instead  of  -|-  Vp  =  0.  In  both  methods  there  is,  in  essence,  a  dual 

velocity  grid  whose  blocks  are  associated  with  the  corners  of  pressure  blocks.  In  contrast, 
as  Sections  2  and  3  will  show,  our  velocity  grid  elements  (control  volumes)  are  associated 
with  the  edges  (or  faces  in  three  dimensions)  of  pressure  blocks,  as  are  the  velocity  degrees 
of  freedom  (e.g.,  Figure  2).  We  also  retain  the  integration  of  k~^  as  in  Eq.  (5),  which  will 
generalize  the  usual  harmonic  averaging  of  in  a  simple  way. 

An  expanded  mixed  method,  reducible  to  a  finite- difference  scheme  by  low-order  inte¬ 
gration,  has  been  formulated,  analyzed,  and  tested  by  Arbogast  ei  al  [3,  4].  The  method  is 
expanded  in  the  sense  that  it  introduces  an  additional  variable  corresponding  to  Vp,  sub¬ 
sequently  eliminating  it  under  some  circumstances.  A  key  to  the  method  is  the  assumption 
that  there  is  a  global  mapping  from  the  irregular  grid  to  a  regular  grid,  which  is  not  the 
case  for  a  grid  such  as  that  in  Figure  1.  At  non-smooth  grid  interfaces,  or  at  interfaces  with 
coefficient  discontinuities,  the  method  must  introduce  Lagrange  multipliers  that  correspond 
physically  to  pressures  on  block  edges  (faces  in  three  dimensions).  With  these  Lagrange 
multipliers,  the  method  obtains  numerically  the  theoretical  convergence  rate  of  for  ve¬ 
locities  at  midpoints  of  edges  [3],  where  h  is  the  grid  size.  The  method  and  its  theory  are 
based  on  the  framework  of  standard  mixed  methods.  Our  method,  based  on  the  alternative 
control- volume  framework,  is  considerably  simpler,  as  its  degrees  of  freedom  are  merely  the 
block  pressures  and  the  edge  fluxes,  with  no  Lagrange  multipliers.  Edge  pressure  values 
(analogous  to  Lagrange  multipliers)  do  appear  in  the  derivation  in  Section  3.2,  but  they  are 
eliminated  before  the  final  system  of  equations  is  reached.  The  resulting  method  shows  a 
convergence  rate  of  for  velocities  (fluxes  across  edges)  in  all  of  the  tests  performed  so  far 
(e.g.,  those  in  Section  4),  except  where  the  solution  is  singular. 

Control  volume  finite  elements.  To  obtain  a  local  discrete  Darcy  law  and  to  avoid 
the  complexities  that  appear  to  be  necessary  with  standard  mixed  methods,  we  consider 
procedures  based  on  control- volume  finite-element  methods  [5,  17].  Such  schemes  have  been 
considered  in  the  petroleum  literature  [14,  20,  13],  but  only  in  a  ‘"point-centered”  framework. 
This  means  that  mass  conservation  is  not  enforced  on  the  blocks  designated  by  the  user,  but 
rather  on  dual  blocks  centered  about  the  vertices  of  the  user  blocks.  Since  the  vertices  are 
presumably  often  located  at  geological  interfaces,  and  since  the  designated  blocks  may  be 
of  significance  to  the  user,  the  most  desirable  approach  would  be  one  that  conserves  mass 
on  the  user  blocks. 

We  briefly  summarize  the  point-centered  approach.  If  Eq.  (1)  is  integrated  over  a  volume 
(area  in  two  dimensions)  V  and  the  Gauss  divergence  theorem  is  applied,  we  obtain 

“  /  —Vp’Xids  =  f  qdx,  (7) 

Jdv  ^  Jv 

where  d\'  is  the  boundary  of  V  and  ds  is  the  measure  on  dV .  A  mesh  of  triangles  is  defined, 
where  p  will  be  computed  at  the  vertices.  With  each  vertex,  one  associates  a  control  volume, 
usually  found  by  taking  the  Voronoi  volume  bounded  by  the  perpendicular  bisectors  of  the 
sides  of  the  triangles.  Then  Eq.  (7)  is  posed  for  each  control  volume,  where  p  is  linearly 
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i+l/2,j 


Figure  3:  Unknowns  and  control  volumes  for  rectangular  grid 


interpolated  from  the  vertices  to  the  interior  of  each  triangle.  Thus,  p  is  represented  by  a 
standard  finite-element  shape  function,  but  instead  of  integrating  against  the  usual  weighting 
functions  one  integrates  over  control  volumes.  This  is  equivalent  to  integrating  against 
weighting  functions  that  are  1  on  one  control  volume  and  0  on  the  others.  Convergence 
theories  [6,  8]  exist  for  quite  general  triangulations. 


2  Rectangular  Grid 

We  formulate  a  control- volume  mixed  finite-element  procedure  for  the  system  Eqs.  (3)-(4), 
using  appropriate  shape  functions  to  represent  the  solution  and  integrating  over  appropri¬ 
ate  control  volumes.  We  first  illustrate  this  for  a  rectangular  grid,  where  the  details  are 
straightforward  and  we  can  use  Raviart-Thomas  elements  as  described  above. 

In  Figure  3,  we  show  typical  unknowns  and  control  volumes.  Let 

—  (^z  — 1/2 5  ^i+1/2)  ^  (yj  — 1/2j  2/;+1/2)5 

Qz-t-l/2j  —  (^zj^z+l)  ^  (^j  — 1/2j  2/j+l/2)5 

—  (^z  — l/2j  ^*+1/2)  ^  + 

where  {xi,yj)  is  the  center  of  the  block  Qij.  As  in  the  standard  mixed  method,  we  associate 
pressure  unknowns  pij  with  block  centers  {xi,yj),  and  a  flux  unknown  (normal  component 
of  V  times  cross-sectional  area  or  length)  with  each  face  (edge  in  two  dimensions).  On  a 
vertical  edge  centered  at  (^*+1/2,  %)?  the  normal  component  is  the  rr-component,  so  we  can 
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denote  the  unknown  by  (/r)z+i/2j-  Similarly,  on  a  horizontal  edge  centered  at  {xi^yj^if2), 
we  associate  The  natural  control  volumes  corresponding  to  (/r)f+i/2j, 

and  {fy)iJ^\|2^  respectively,  are  Qi  j^  Qi^ii2,3,  and  Qij^ii2,  as  seen  in  Figure  3.  From  a 
physical  point  of  view,  we  can  think  of  Qf+i/2j  as  a  “tank”  with  pressures  pij  and 
imposed  at  the  two  ends,  and  similarly  for  Qij^if2> 

We  wish  to  integrate  Eq.  (3)  over  control  volumes  to  obtain  local  discrete  Darcy  laws 
on  the  “tanks.”  Note  that  Eq.  (3)  is  a  vector  equation  that  can  be  resolved  into  the  two 
components 


dp 

(8) 

Vy  +  ~  =  0, 

dy 

(9) 

where  A  —  k/ fi.  We  integrate  Eq.  (8)  and  Eq.  (9)  over  and  Qij-j.i/2,  respectively, 

because  these  equations  involve  Vx  and  Vy.  Integrating  out  the  partial  derivatives  of  p,  this 
yields 


rv 

J  Xi  J  y 
fyj+i 

Jyj  Jxi 


yj+i/2 


A  ^Vx{Xjy)  dy  dx + 


^yj-i/2 

pj  +  l  f^t  +  l/2 
■1/2 


/ 


yj+i/2 


ip(xi+i,y)  -p(xi,y))dy 


yj-i/2 


A  ^Vy(x,y) 


dxdy+  /  {p{x,yj+ 


i)  -  p(x,yj))dx 


0,  (10) 

0.  (11) 


Using  the  shape  functions  forp,  and  Vy  as  described  above,  the  integrals  in  Eqs.  (lO)-(ll) 
are  readily  expressed  in  terms  of  the  unknowns  pij,  pi+ij,  pij+i,  ifr)i-i/2,j^  {fx)i+i/2,j, 
{fx}i+3/2,j ,  {fy)i, }-i/2,  ify)i,j+i/2y  and  {fy)ij+z/2-  We  then  obtain  the  discrete  Darcy 
equations  on  the  “tanks”:  in  the  a;-direction  on  Qi+i/oj, 


where 


<^i+l/2,j-,i-l/2,jifx)i-l/2J  +  0'i+l/2,j;i+l/2,j(fx)i+l/2,j 

+(‘i+l/2J;i+3/2,j{fx)i+3/2,j  +Pi+lJ  -  Pi,j  =  0,  (12) 


1  \,j  ■> 

h  +  l/2j;!-l/2j  = 

(13) 

h+l/2,j;i+l/2j  = 

8  IQi.il  ~  *^*-1/2) 

+ 

8  "^i+l.i  /  ^2 

(14) 

■i+l/2,j-i+3/2,j  = 

1  A~1 

^  /  n2 

(15) 

and  in  the  j/-direction  on  Qij+1/2, 

(^i,j  +  l/2;i,j-l/2{fy)i,j-l/2  +  ^i,j  +  l/2-,i,j  +  l/2ify)iJ  +  l/2 

+^i,j+l/2-,i,j+3/2{fy)i,j+3/2+Pi,j  +  l  -  PiJ  =  0,  (16) 
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where  coefficients  in  Eq.  (16)  are  obtained  by  equations  analogous  to  Eqs.  (13)-(15).  On 
each  block  Qij,  p  is  constant,  /j,  varies  linearly  with  x  and  is  constant  in  y,  and  fy  is 
constant  in  x  and  varies  linearly  with  y. 

We  also  integrate  Eq.  (4),  this  time  over  the  control  volumes  Qij.  Applying  the  Gauss 
divergence  theorem  to  convert  the  left-hand  side  into  a  boundary  integral  (four  edge  inte¬ 
grals),  we  have 


fyj+i/2  yJ/j+i/i  r®i+i/2 

/  Vr{xi+i/2,y)dy  -  /  v^{xi_i/2,y)dy+  /  Vy{x,yj+i/2) 

Jyj-1/2  ^yj-if2  ^^1-1/2 

r<^t+i/2  ryj+i/2  r^i+i/2 

-  vy{x,yj^i/2)dx=  /  q{x,y)dxdy. 


dx 


(17) 


Again,  the  integrals  are  expressed  in  terms  of  {fy)i,j-i/2y  and 

ify)ij+i/2i  yielding  the  discrete  mass  conservation: 

{fx)i^l/2,j  -  {fx)i+l/2,j  +  {ftj)i,j-l/2-  {fy)i,j  +  l/2  —  “IQijkij*  (1^) 


Eqs.  (12),  (16),  and  (18)  thus  give  rise  to  a  symmetric  system  of  linear  equations  that  is 
solved  for  the  pressures  at  block  centers  and  the  fluxes  across  edges: 


■  0  ■ 

fx 

0 

0  My  Ny 

fy 

0 

.  Nj  Nj  0  _ 

P 

where  M^-  is  tridiagonal  (for  unknown  ordering  by  horizontal  rows,  so  that  (i  -  1/2,  j), 
(i+  1/2,  j),  (i  +  3/2,/)  are  consecutive),  My  is  tridiagonal  (ordering  by  vertical  columns 
so  that  {ij  -  1/2),  (ij  +  1/2),  {ij  +  3/2)  are  consecutive)  and  each  row  of  and  Ny 
contains  one  +1  and  one  -1,  corresponding  to  the  two  adjacent  block  pressures. 

It  is  instructive  to  relate  this  control- volume  mixed  finite-element  method  to  the  familiar 
block-centered  finite-difference  approach.  Eq.  (17)  would  be  the  usual  block-centered  mass- 
balance  equation  if  the  normal  velocities  on  edges  were  given  by  discrete  pressure  gradients 
multiplied  by  harmonically  averaged  mobilities.  Examining  Eqs.  (lO)-(ll),  we  see  that  this 
would  be  the  case  if  the  and  Vy  integrals  approximated  and  Vy  by  constants  on  their 
respective  control  volumes  (or,  equivalently,  if  a  midpoint  integration  rule  were  used).  In 
matrix  terms,  M.r  and  would  become  diagonal,  so  the  system  would  have  the  form 


■  M' 

N  ■ 

■  /  ■ 

0 

0 

.  P  . 

.  -\Q\<i . 

with  M'  diagonal,  and  elimination  of  /  would  yield 

N^M'“^Np=  \Q\q. 


(20) 


(21) 


This  demonstrates  the  close  relationship  between  the  control-volume  mixed  method  and 
block-centered  finite  differences.  Both  methods  involve  local  mass  conservation  and  a  local 
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(0.0) 


(1.0) 


(X  .y  ) 
11  11 


Figure  4:  Reference  quadrilateral  Q  and  quadrilateral  Q 


Darcy  law.  The  mixed  method’s  higher-order  approximation  of  and  Vy,  each  varying 
linearly  in  its  own  direction,  couples  the  velocities  in  Eqs.  (10)”(11)  and  makes  the  methods 
different.  This  increases  the  accuracy  of  the  solution  (for  example,  see  Problem  1  of  Section 
4),  but  makes  it  more  expensive  to  solve  the  discrete  system. 

It  is  also  useful,  for  the  sake  of  generalization  to  irregular  grids,  to  take  a  vector  point  of 
view  of  Eqs.  (lO)-(ll).  In  passing  from  Eq.  (3)  to  Eq.  (8),  we  took  the  ^‘-component  of  Eq. 
(3);  this  is  equivalent  to  taking  the  dot  product  of  Eq.  (3)  with  the  unit  vector  x  =  (1,  0).  In 
Eq.  (10)  we  restricted  the  integration  to  the  control  volume  (5z+i/2,j-  Thus,  we  can  obtain 
Eq.  (10)  by  taking  the  dot  product  of  Eq.  (3)  with  a  vector  field  that  is  x  on  Qi-\-i/2,j  and 
zero  elsewhere,  then  integrating  over  Q.  This  vector  field  is  the  finite-element  vector  “test 
function”  corresponding  to  Eq.  (10).  Similarly,  Eq.  (11)  relates  to  a  vector  test  function 
that  is  y  =  (0, 1)  on  and  zero  elsewhere.  This  perspective  seems  pointless  on  a 

rectangular  grid,  where  components  are  easy  to  work  with,  but  it  will  be  helpful  otherwise. 


3  Irregular  Quadrilateral  Grid 

With  the  rectangular  case  as  a  guide,  we  develop  a  formulation  for  general  quadrilaterals. 
One  important  step  is  to  be  able  to  relate  a  general  quadrilateral  to  a  reference  one.  Consider 
the  quadrilateral  Q  in  Figure  4,  which  is  assumed  to  have  vertices  at  (^oOi2/oo)j  (^oi,2/oi), 
(^io,yio),  and  (xii,t/ii).  Let  the  reference  quadrilateral  Q  be  the  unit  square.  There  is  a 


Figure  5:  Velocity  basis  function  on  quadrilaterals 


unique  bilinear  mapping  of  Q  onto  Q  that  sets  up  coordinates  on  Q: 


x{x,y)  —  soo  +  (a;io  —  a;oo)^  +  (a^oi  -  a-’oo)y 

+(^11  -  a:io  -  a-'oi  +  xoo)xy, 

(22) 

2/(i ,  y)  -  yoo  +  (yio  -  j/oo)*  +  (j/oi  -  yoo)y 

+(yii  -  yio  -  yoi  +  yoo)xy. 

(23) 

The  resulting  coordinate  lines  are  as  depicted  in  Fig.  4.  The  pressure  in  Q  is  associated 
with  the  image  of  the  center  of  Q,  i.e.,  with  the  node  ?/(|,  ^))  indicated  by  x  in 

the  figure.  Note  that  this  is  not  generally  the  centroid  of  Q.  As  long  as  Q  is  a  convex 
quadrilateral  (all  angles  less  than  180  degrees),  the  bilinear  mapping  has  an  inverse.  We 
assume  henceforth  that  the  quadrilaterals  are  convex  so  that  for  each  (x’,  y)  £Q  the  inverse 
mapping  gives  an  associated  G  Q .  There  is  thus  a  one-to-one  correspondence  between 

points  in  the  physical  space  Q  and  the  reference  space  Q. 

3.1  Shape  functions  and  unknowns 

Now  consider  the  extension  of  the  control- volume  mixed  formulation  to  general  quadrilater¬ 
als.  To  maintain  continuity  of  flux,  we  want  the  normal  component  of  a  velocity  function  to 
be  constant  on  each  edge.  Then  we  can  associate  degrees  of  freedom  with  fluxes  on  edges, 
as  in  the  rectangle  and  triangle  cases.  In  Figure  5  we  show  two  adjacent  quadrilaterals  with 
the  coordinates  determined  by  the  mapping  Eqs.  (22)-(23).  The  velocity  vector  function 
that  has  normal  component  1/|F'|  on  the  common  edge  of  length  \E\  (hence  has  flux  1) 
and  0  on  the  other  edges  is  pictured.  It  is  oriented  along,  say,  ar-coordinate  lines,  and  has 
constant  normal  component  on  each  complementary  y-line,  with  the  magnitude  of  the  flux 
varying  linearly  in  the  a?-direction.  We  now  describe  this  vector  function  analytically. 

First  we  identify  significant  directions  in  the  quadrilateral.  Referring  to  Eqs.  (22)-(23), 


define 


X(£,p)  = 


Y(x,y)  = 


\5i;  ’  d£  J 

(xio  -  a;oo  +  (a^ii  -  a^io  -  a;oi  +  a;oo)y, 
2/10  -  2/00  +  (2/11  -  2/10  -  2/01  +  2/oo)y), 

[dy’dyj 

(xoi  -  *00  +  (a;ii  -  xio  -  a;oi  +  xoo)£, 
2/01  -  2/00  +  (2/11  -  2/10  -  2/01  +  2/oo)£)- 


(24) 


(25) 


These  can  be  viewed  as  the  images  of  the  vectors  (1,0)  and  (0, 1),  respectively,  under  the 
mapping  from  Q  to  Q.  We  have  defined  them  for  (i,  y)  in  the  reference  quadrilateral,  but 
because  of  the  one-to-one  correspondence  mentioned  previously,  we  can  just  as  well  consider 
them  defined  for  (x,  y)  in  the  physical  quadrilateral.  In  the  physical  space,  they  point  in  the 
directions  of  the  coordinate  lines  pictured  in  Figure  4.  However,  they  are  not  unit  vectors; 
their  length  depends  on  the  size  of  Q  and  they  have  the  dimensions  of  length.  Define  also 
the  corresponding  unit  vectors  and  normal  vectors; 


Here  x  and  y  are  unit  vectors  in  the  directions  of  X  and  Y,  respectively,  n^,  is  a  unit  vector 
normal  to  Y,  and  iij,  is  a  unit  vector  normal  to  X.  Figure  5  shows  x  and  y,  while  Figure  6 
shows  nj,  and  iiy . 

Returning  to  the  vector  function  v  in  Figure  5,  let  Q  be  the  left-hand  quadrilateral.  To 
evaluate  v  at  {x,y),  first  use  the  inverse  mapping  to  find  the  corresponding  {£,y).  Then 
v(x,y)  is  the  vector  in  the  direction  of  X  whose  nj--component  (i.e.,  v-iij,)  is  equal  to  x/||Y|| 
(so  that  the  flux  across  the  “vertical’’  line  through  (x,y)  is  £).  After  some  manipulation, 
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Figure  6:  Control- volume  mixed  finite  elements  on  quadrilaterals 


one  can  show  that  this  is  given  by 

v{x,y) 


xX 

J(x,y)’ 


where 

dx  dy  dx  dy 

—  Qx  dif  dy  dx 

-  [(a^io  “  a?oo)(yoi  -  yoo)  -  (^01  -  a^oo)(yio  -  Voo)] 

+[(^’10  -  xoo){yn  -  yoi)  ~  (^ii  -  i^oi)(yio  -  yoo)]^ 
-{-[{xii  -  xio){yoi  -  yoo)  -  (^01  -  a^oo)(yii  -  yio)]y 


(30) 


(31) 


is  the  Jacobian  of  the  mapping  in  Eqs.  (22)-(23)  from  Q  to  Q.  Note  that  J{x,y)  = 
||X||||Y||sin6f,  where  0  is  the  angle  between  X  and  Y.  For  a  rectangular  grid,  this  an¬ 
gle  is  90  degrees,  so  that  n^:  —  x  and  we  obtain  v(j^,  y)  ===  :cx/||Y||,  as  we  would  expect.  In 
the  right-hand  quadrilateral,  everything  is  the  same  except  that  1  —  x  replaces  x.  Similar 
expressions  hold  for  a  velocity  basis  function  corresponding  to  a  “horizontal”  edge,  with  x, 
X,  and  X  replaced  by  y,  Y,  and  y,  respectively. 

The  unknowns  that  we  use  to  describe  the  velocity  are  the  fluxes  across  edges,  namely 
\Ei+i/2,j\{v''^x)i+i/2j  and  j+i/2|(v  •  ny)£j+i/2,  which  we  abbreviate  to  {fx)i+i/2j  and 
{fy)ij^i/2  analogy  with  the  rectangular  case.  These  velocity  shape  functions  and  un¬ 
knowns  can  be  obtained  from  those  on  rectangles  by  a  so-called  Piola  transformation  [7]. 
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For  pressure,  the  natural  choice  of  shape  functions  is  still  piecewise  constants,  introducing 
no  additional  complications,  and  the  unknowns  are  pij  as  in  the  rectangular  case. 


3,2  Test  functions  and  control  volumes 


To  obtain  discrete  equations  from  which  we  can  solve  for  pressures  at  block  centers  and  fluxes 
across  edges,  we  must  choose  suitable  control  volumes  and  mimic  the  integrations  leading 
to  Eqs.  (10)“(11)  and  (17).  For  the  integrations  of  Eq.  (3),  we  use  images  of  rectangular 
control  volumes  Qi^i/2,j  and  Q*  j-(-i/2  under  the  mapping  in  Eqs.  (22)-“(23).  This  is  pictured 
in  Figure  6,  where  the  control  volume  associated  with  the  common  edge  consists  of  the 
image  of  (^,  1)  x  (0, 1)  under  the  mapping  to  the  left-hand  quadrilateral  and  the  image  of 
(0,:|)  X  (0,1)  under  the  mapping  to  the  right-hand  quadrilateral.  In  physical  space,  this 
can  be  described  by  taking  the  midpoints  of  the  four  edges  adjacent  to  the  common  edge  in 
question,  then  joining  the  two  pairs  of  midpoints  by  line  segments.  We  denote  such  control 
volumes  by  (5z+i/2j  and  as  we  did  previously  for  the  rectangular  grid.  Qi+i/nj 

in  Figure  6  will  be  the  “tank”  with  pressures  pij  and  Pi^ij  at  the  two  ends.  For  the 
integrations  of  Eq.  (4),  we  simply  take  the  quadrilateral  blocks  Qij  as  control  volumes. 

Continuity  equation.  For  the  integrations,  we  also  require  test  functions.  For  Eq.  (4), 
these  are  simply  scalar  characteristic  functions  of  the  control  volumes,  i.e.,  functions  that 
are  1  on  one  volume  and  zero  elsewhere.  If  we  denote  the  edges  of  Qj  j  by  etc., 

and  integrate  Eq.  (4)  over  the  control  volume,  the  Gauss  divergence  theorem  (using  the  fact 
that  the  normal  velocity  component  is  constant  on  each  edge)  yields 


(v  •  n^)i+i/2j|^i+i/2,j|  -  (v  ’  nj:)j_i/2j\Ei_i/2,j\  +  (v  ^  ny)ij^i/2\Eij^i/2\ 

—  (v  •  n^ j_i/2|  =  /  qdzj  (^*^) 

so  that  Eq.  (18)  is  obtained  for  quadrilaterals  as  well  as  rectangles,  where  dz  is  the  two- 
dimensional  measure  on  Qij^  The  equation  is  easily  incorporated  into  the  discrete  system 
as  before. 

Darcy  equation.  For  Eq.  (3),  the  situation  is  more  complicated.  On  the  rectangular 
grid,  the  test  function  given  by  x  on  Qi^ifoj  and  zero  elsewhere  allowed  the  .r-partial 
derivative  of  77  to  be  integrated  out,  leaving  integrals  of  p  on  lines  that  were  interior  to  the 
constant-pressure  blocks.  This  is  the  desired  outcome,  and  we  show  how  to  achieve  this  on 
a  general  quadrilateral  grid. 

Let  Qi^i/4j  and  Qi4.sf4j  denote  the  “left-hand  half”  and  “right-hand  half,”  respectively, 
of  Qi^i/2,j-  Then  (5?:+i/4j  is  the  image  of  the  right-hand  half,  (1/2, 1)  x  (0, 1),  of  Q  under 
the  mapping  Eqs.  (22)-(23).  Using  X  as  the  test  function,  the  p  integral  analogous  to  the 
one  in  Eq.  (10)  is 


[  Vp-Xdz  =  f 

-  L 


+1/4, > 


dpdx  dpdy^. 
dx  dx  dy  dx  ^ 
dp 

-^dz 

ox 
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(33) 


where  J  =  Jij  is  the  Jacobian  of  the  mapping  from  Q  to  as  in  Eq.  (31).  Since  J  is 
linear  in  f ,  dJjdx  =  6  is  constant,  and  Eq.  (33)  becomes 

/  VpXrfz  =  Aj(l.y)p(l,y)-J(l/2,y)p(l/2,y)]Jy 

“Si 

Now  we  suppose  that  p  is  approximated  by  a  linear  (not  bilinear)  polynomial  in  x  and 
y,  an  approximation  of  error  O(h^),  where  h  is  the  diameter  of  Qi^if4j.  This  will  allow  us 
to  reduce  the  above  expression  to  an  appropriate  numerical  scheme.  Since  J  is  also  linear, 
the  Jp  integrals  can  be  evaluated  by  Simpson’s  rule  in  y: 

Jil,y)pil,y)dy  =  ^Jil,0)pil,0)  +  pil,l/2)pil,l/2) 

J(l/2,  y)p(l/2,  y)  dy  =  ^(1/2,  0)p(l/2, 0)  +  ^(1/2,  l/2)p(l/2, 1/2) 

+  A(1/2,1)K1A1)-  (35) 

For  the  p  integral,  we  use  the  trapezoidal  rule  in  x  and  Simpson’s  rule  (higher  order  than 
necessary,  but  easier  to  combine  with  the  other  terms)  in  y: 

Si  "  s|f«l/2,0)  +  p(l,0))  +  lg(,,(l/2,l/2)  +  p(l,l/2)) 

+  ^|f(p(l/2,l)  +  p(l.l))-  (36) 

Substituting  Eqs.  (35)-(36)  into  Eq.  (34)  and  collecting  coefficients,  we  have 

-(jW2,0)+1|)p(1/2,0) 
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-(|j(l/2,l/2)+ig)p(l/2,l/2) 

=  ^J(3/4,0)(p(1.0)-p(l/2,0)) 

+1  J(3/4,  l/2)(p(l,  1/2)  -  p(l/2, 1/2)) 
+^J(3/4,l)(p(l,l)-p(l/2,l)).  (37) 

where  the  last  step  uses  the  linearity  of  J  with  respect  to  x.  With  p  being  linear,  dp/dx  is 
constant,  so  that  the  three  p  differences  just  obtained  are  all  equal;  hence,  using  the  linearity 
of  J  with  respect  to  y,  Eq.  (37)  reduces  to 


^*  +  1/4,/ 


Vp-Xdz-  Jij(3/4,  l/2)(pi+i/2j  -  Pi,j), 


(38) 


recalling  that  pij  =  p{l/2, 1/2)  and  letting  Pi+i/oj  =  p{l,  1/2)  be  a  pressure  value  on  the 
edge  The  value  Pi^i/2j  is  not  one  of  the  desired  block-center  pressure  unknowns 

and  we  wish  to  eliminate  it.  A  similar  derivation  for  leads  to 


(39) 


Hence,  by  choosing  the  test  vector  field 

[  X/Aj(3/4, 1/2)  on(5,+i/4j, 

y^i+i/2j  =  <  X/J^+i,,-(l/4,l/2)  0TiQi^s/4,j,  (40) 

[  0  elsewhere, 

we  combine  constant  multiples  of  Eqs.  (38)-(39)  into 

/  Vp-  w,+i/2j  dz  =  (41) 

and  the  edge  value  p*+i/2j  has  been  eliminated.  Note  that  we  did  not  have  to  require  that 
p  be  piecewise  constant  in  this  derivation,  though  we  will  generally  think  of  the  numerical 
approximation  of  p  in  this  way. 

The  step  just  completed  is  the  elimination  of  the  analogues  of  the  Lagrange  multipliers 
of  Arbogast  ei  aL  [3,  4],  mentioned  in  the  discussion  of  other  methods  for  irregular  grids 
in  Section  1.  The  ability  to  carry  out  this  step  is  a  special  property  of  the  control- volume 
mixed  method,  as  opposed  to  the  standard  framework  in  which  the  vector  shape  and  test 


14 


functions  are  the  same.  In  the  latter  case,  the  test  function  must  have  continuous  normal 
flux,  and  there  is  no  freedom  to  choose  weights  as  in  Eq.  (40)  above.  There  is  no  physical 
reason  for  this  constraint  (unlike  flux  continuity  for  the  shape  functions),  which  is  an  artifact 
of  the  numerical  approach.  The  control- volume  formulation  allows  (indeed,  compels)  flux 
discontinuities  in  the  test  functions. 

At  this  point  we  have  chosen  a  test  function  and  have  integrated  the  p  term  of  Eq.  (3). 
For  the  v  term,  we  consider 

/  A“^v-w<+i/2,jrf2  =  /  •X/Jij(3/4,  l/2)(/z 

JQ.  + 

+  /  A-^\  .vXM+i,,(l/4,l/2)rfz,  (42) 

where  A  may  be  a  full  anisotropic  tensor.  To  find  the  desired  coefficients,  write 

V  =  (/ar)*--l/2,jVj_l/2j  +  +  (/a7)2-{-3/2jVj+3/2,j 

+  (/2/)2,j  +  l/2Vzj  +  l/2  +  (/y)*j-l/2Vi,j_l/2 

+  (/y)z+l,j  +  l/2V*  +  ij  +  i/2  +  (/^)i+l,j-l/2V*  +  lj-l/2 

+  other  terms,  (43) 

where  (for  example)  Vi+i/2j-  is  the  velocity  field  with  flux  1  across  F^2+i/2,j  and  0  across  all 
other  edges.  Then,  substituting  for  v,  we  obtain  the  discrete  Darcy  equation  analogous  to 
Eq.  (12) 


«2-(-l/2j;«-l/2,j(/rr)/-l/2,i  d"  «*+l/2,i ;z  +  l/2,j (/rr )2+l/2,i  d"  ^2+l/2j;i+3/2j(/r)i+3/2,i 
d-a2  +  i/2j;*,j+l/2(/y)zj  +  l/2  +  «i  +  l/2j;z,j~l/2(/y)ij-l/2 
+af+l/2j;2  +  lj+l/2(/y)2  +  l,j  +  l/2  +  «?+l/2,i;i+l  j -l/2(/2/ )z+lj-l/2 

-Pij  =  0,  (44) 

where  (for  example  in  analogy  with  Eq.  (14)) 

Cti+l/2,j;i  +  l/2j  =  /  ^2  Jvi+i/2j  ■  X/Jij(3/4,  1/2)  dz 

dQ^+l/4,j 

+  /  A“,\^  •  X/ J2+i.,-(1/4,  1/2)  dz.  (45) 

To  evaluate  the  integral  in  Eq.  (45),  first  note  that  is  parallel  to  X.  Thus, 

the  unit  vectors  V2+i/2j7l|vi4-i/2j||  ^^d  X/||X||  point  in  the  same  direction,  and  are  there¬ 
fore  equal.  Also  note  that  J  —  j|X||||Y||(vj_|_i/2j  •  '^x)/\\^i~\-\f2,3\\^  because  the  angle  0 
between  X  and  Y  (hence  between  V2+i/2j  and  Y)  is  the  complement  of  the  angle  7]  be¬ 
tween  Vj_|.i/2y  and  so  that  sin^  —  cost?  =  (v2+i/2j  *  Jajp)/||vj4.i/2j||.  Furthermore, 
l|Y|l(  Vj_|_i/2j  •  nj;)  is  the  flux  that  varies  linearly  across  (5i+i/4j,  being  equal  to  1  at  the 
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right  edge  and  1/2  at  the  left  (the  ‘Vertical”  center  line  of  Qij);  hence,  the  flux  is  equal  to 
X.  Finally,  recall  that  Av  •  w  =  v  •  A^w  if  A  is  a  square  matrix  with  transpose  A^.  Thus: 


■  ^dz 


'Q,+i/4.j  l|V,  +  l/2j|| 


'Qt-+l/4,j 


l|X||||Y||(v  i  +  l/2,j  ■  Px)  X  /.-IxTy^ 

j  ||xn 


/  ||Y||(vi+i/2,;-n,)(A-/X)-xV.' 

I  I  x(K}^)-^dxdy, 

Jo  Jl/2 


where  the  last  step  is  a  change  of  variable.  Similarly,  we  find  that 


'Qt  +  3f4,j 


rl  .1/2 

{l-i)(Ar^,jX)-Xdxdy. 


Hence,  combining  Eqs.  (45)-(47), 


=  ^,7574-/2)/  I  i(A-‘X)^Xdi:dii 


1  pl/2 


<^i+ij(l/4i  1/2)  Jo  Jo 


il-x)iA-^,jX)-Xdxdy.  (48) 


Next,  consider  v^_i/2j  in  order  to  obtain  02+1/2, in  analogy  with  Eq.  (13).  Of  the 
two  “halves”  j  and  (5^+3/4 j  where  does  not  vanish,  Vi-i/nj  is  nonzero  only 

on  Qi^if4  j,  Reasoning  as  above, 


Similarly,  'Vi+s/nj  is  nonzero  only  on  <5,:+3/4j,  and  Eq.  (15)  has  the  analogue 


/  /  (l-i)(A-/X).Xdi-#. 

Jo  Jl/2 


h-+l/2j;i+3/2,j 


,,■(1/4,  1/2)  i  i  i(Ar+\,,.X).Xdird2/. 


The  logic  is  slightly  different  if  the  term  for  a  horizontal  edge,  e.g.  Vj j+1/2  to  obtain 
02+1/2, j;2,/+i/2t  is  considered.  These  terms  were  0  in  the  rectangular  case.  Here  again 
only  one  of  the  “halves”  is  relevant,  in  this  case  (3,:+i/4j.  Now  v^j+i/2  is  parallel  to  Y, 
so  that  v,-.,-+i/2/||vi,, -+1/211  =  Y/||Y||,  and  J  =  ||X||||Y||(v+, -+1/2’ •  ny)/||v+, +1/211;  also, 
||X||(v,,/+i/2  •  iiy)  =  y.  Then  Eq.  (46)  is  replaced  by 


d^i,jX,j+i/2  Xdz 


-  L 


IK,i+l/2||Tp^-(A-/)^Xd3 
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(51) 


=  /  ||X|i(v,,,.,i/2-%)(A-/Y)-xlrfz 

=  t  t  y{kT]Y)^lLdxdy, 

Jo  Jl/2 

SO  that  ^  ^ 

ljiA7]Y)-Xiidy^ 

By  analogous  steps,  the  other  three  coefficients  are: 

1 

a,+l/2.;;.+l,i+l/2  =  . (1/4^  1/2)  ^  Jo  ^(^r+i.iY)-Xd£rfy, 


1  /•!  /•1/2 
[TJrwii  (») 


(52) 

(53) 

(54) 


ai+l/2,j;t  +  l,i-l/2  -  J.^j^^.(l/4 

This  completes  the  description  of  the  Darcy  equation  for  the  vertical  edge  £'j+i/2,j.  The 
Darcy  equation  for  the  horizontal  edge  Eij^i/2, 


ai,i  +  l/2;i,j-l/2(/y)»j-l/2  +  ai,i  +  l/2;i  j  +  l/2(/y  )i,j  +  l/2  +  ai.j  +  l/2;ij+3/2(/y)i,i+3/2 
+  0,i,j  +  l/2-,i  +  l/2,jifx)i+l/2,j  +  (^i,j  +  l/2-,i-l/2,j{fx)i-l/2,j 
+  ai,j  +  l/2;>  +  l/2j  +  l(/x)i+l/2,j+l  +  «!,j+l/2;j-l/2,i  +  l(/®)!-l/2,l+l 

+Pij  +  1-Pi,j  =0,  (56) 


is  derived  in  a  completely  analogous  fashion.  The  coefficients  are  defined  by  equations 
similar  to  Eqs.  (48)-(50)  and  (52)-(55). 

Assuming  that  the  reciprocal  mobility  A" Ms  a  constant  tensor  on  each  grid  block,  the 
integrals  in  the  a  coefficients  are  straightforward  to  evaluate  analytically.  The  dot  products 
in  these  integrals  are  simply  quadratic  polynomials  in  x  and  y  (total  degree  2,  so  that  the 
highest-order  terms  are  xy,  y^).  Explicit  expressions  for  the  a-coefficients  can  be  found 
in  [21].  These  can  be  evaluated  once  and  stored  for  use  throughout  the  life  of  the  grid 
block  in  the  simulation.  Even  in  multiphase  or  variable- viscosity  flow,  where  A”^  is  time- 
dependent,  the  dependence  is  restricted  to  a  scalar  multiple  of  the  tensor,  so  that  the  above 
double  integrals  can  be  stored  and  later  multiplied  by  a  variable  scalar. 


3,3  Discrete  system  of  linear  equations 

The  discrete  system  of  the  control- volume  mixed  finite  element  method,  with  irregular 
quadrilateral  grid  and  full  anisotropic  tensor  permeability,  consists  of  ^Vertical” -edge  Darcy 
Eq.  (44),  “horizontal” -edge  Darcy  Eq.  (56),  and  continuity  Eq.  (18).  The  non-symmetric 
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linear  system  can  be  written  in  the  form 


IVIyj: 

yiyy 

N, 

l-iy 

0 

fx 

fy 

P 


(57) 


where  and  are  tridiagonal  with  the  same  nonzero  pattern  as  in  the  rectangular 
case  Eq.  (19)  (though  with  different  values),  and  and  Ny  have  the  same  ±1  entries  as  in 
Eq.  (19).  Mj-y  and  each  have  four  nonzero  bands,  corresponding  to  the  four  nonzero 
a-coefBcients  involving  (A“^  Y)  •  X  for  and  (A-iX).Yfor 

Note  that  the  a-coefficients  in  Eqs.  (48)“(50),  (52)-(55)  amalgamate  the  complexities  of 
distorted  grids  and  of  anisotropic  tensor  permeabilities.  With  a  scalar  permeability  and  an 
orthogonal  grid,  one  sees  that  a  dot  product  such  as  (A~^X)  •  Y  vanishes,  because  A‘^X 
is  parallel  to  X  (due  to  the  scalar  A”^)  and  hence  perpendicular  to  Y  (by  orthogonality). 
Then  the  matrix  M  of  a-coefficients  reduces  to  tridiagonal  form,  as  observed  previously 
for  the  rectangular  case.  If  either  condition  fails  to  hold,  (A"^X)  ■  Y  can  be  nonzero,  and 
additional  bands  can  appear  in  M.  The  nonzero  pattern  with  both  distortion  and  anisotropy 
is  the  same  as  with  either  one  alone.  Within  the  bounds  of  consistent  discretization,  the 
expressions  in  Eqs.  (18),  (44),  (48)-(50),  (52)-(56)  seem  as  simple  as  one  could  reasonably 
hope.  Given  a  tensor  A“\  there  is  a  theoretical  possibility  of  choosing  distorted  X  and  Y 
such  that  (A  ^X)  •  Y  vanishes,  resulting  in  a  sparser  M,  but  the  practical  significance  of 
this  is  not  clear. 


3.4  Extension  to  Three  Dimensions 

It  is  important  to  realize  that  the  system  obtained  here  on  quadrilaterals  extends  readily  in 
three  dimensions  to  hexalaterals  H  that  are  trilinear  images  of  a  unit  cube  H  =  [0, 1]^.  The 
faces  of  such  hexalaterals  may  not  lie  in  a  plane,  but  this  is  not  a  concern  in  principle  because 
the  curvilinear  faces,  normal  vectors,  and  fluxes  are  uniquely  determined.  A  “horizontal”- 
face  Darcy  equation  similar  to  Eq.  (44)  or  Eq.  (56)  would  have  the  form 

j,A-+l/2;2,i,A:-l/2(/x  )ij,A*-l/2  T  <^ij,A'+l/2;i  j,^-)-3/2(  A  )f  j-,A:-f  3/2 
+  a^J,^'+l/2;^4-l/2J^A;(/a;)^  +  l/2J,A:  +  «i,i,A-+l/2;?:-l/2j,A-(/j*)7;-l/2j,fc 
+«i,i,fc  +  l/2;2  +  l/2J.^+l(/F)i+l/2,j,A*+l  +  J  ,A:  +  l/2;z- 1/2  j,A-  +  l(/r  )*- 1/2  j,A:  +  l 

+^^^J,A♦+l/2;^J  +  l/2,A’(/y)f  J  +  1/2,A;  +  (2, k{fy)i,j -1/2, k 

+  ^z  j,A:  +  l/2;f  j+l/2,A:+l(/y  )z ,i  +  l/2,AT  +  l  +  ^^i,A:+l/2;^  j  - l/2,A;+l(/2/ )?:  j  -  1/2,A:  +  1 

+  1  ^PiJ,k  =  0,  (58) 


with,  for  example, 


^i,i,fc+l/2;2  +  l/2j,A: 


III  HKlk^)‘^dxdydz. 

J1/2J0  Jo 


Jij^ki^/"^^  1/2, 3/4)  J 1/2  Jo 
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The  Darcy  equations  for  “vertical”  faces  normal  to  x-fluxes  and  ^/-fluxes  would  be  obtained 
analogously.  The  continuity  equation  would  have  the  form  of  Eq.  (18)  with  two  additional 
/j  terms,  one  with  each  sign.  We  note  that  the  control-volume  mixed  finite  element  method 
has  been  used  in  three-dimensional  astrophysical  applications  [9].  However,  only  Cartesian 
grids  were  considered  and  the  mobility  coefficient  A  was  constant  and  scalar. 

3.5  Relation  to  Block-Centered  Finite  Differences 

On  rectangles,  we  saw  that  the  control- volume  mixed  method  reduced  to  block-centered  fi¬ 
nite  differences  if  the  normal  fluxes  were  constant  instead  of  linearly  varying.  The  equivalent 
reduction  in  the  present  setting  is  to  replace  the  factors  x,  1  —  x,  y,  or  l  —  y  in  Eqs.  (48)-(50), 
(52)-(55)  by  1  if  they  are  greater  than  1/2,  and  by  0  if  they  are  less  than  1/2.  Then  Mxr  and 
Myj,  become  diagonal  (Eqs.  (49)-(50)  yield  zeros),  but  M^y  and  Myx  do  not  vanish.  This 
reflects  the  necessity  of  retaining  cross-derivative  information  in  a  consistent  approximation 
when  grid  distortion  or  anisotropy  is  present.  The  inconsistency  of  corner-point  geometry 
[18]  is  a  consequence  of  its  suppression  of  this  information,  so  as  to  work  within  a  5-point 
stencil  in  two  dimensions  and  a  7-point  stencil  in  three.  If  distortion  and  anisotropy  are 
mild,  then  M  is  strongly  diagonally  dominant,  and  it  should  be  reasonable  to  approximate 
M~i  by  a  matrix  with  the  same  sparsity  pattern  as  M.  The  resulting  analogue  of  Eq.  (21) 
would  have  a  9-point  stencil  in  two  dimensions  and  a  19-point  stencil  in  three,  the  same 
connection  structure  found  by  Arbogast  ei  al  [4].  This  has  not  been  implemented  at  this 
writing  and  will  not  be  discussed  further  here. 


4  Results 


The  control-volume  mixed  finite  element  method  has  been  tested  on  a  variety  of  two- 
dimensional  problems,  involving  uniform  and  irregular  grids,  scalar  and  tensor  permeabil¬ 
ities,  and  constant  and  variable  permeabilities.  The  velocities  exhibit  second-order  con¬ 
vergence  in  all  situations  except  where  the  exact  solution  has  a  singularity,  in  which  case 
second-order  convergence  is  not  possible.  Following  is  a  representative  sampling  of  these 
results.  Additional  results  can  be  fonnd  in  [21]  and  [16]. 

In  most  of  the  test  problems,  an  analytical  solution  was  known.  Otherwise,  a  suitable 
fine-grid  numerical  solution  was  used  for  this  purpose.  Let  p  and  (u®,  Vy)  denote  the  exact 
pressure  and  velocity,  respectively,  with  P  and  (14,14)  being  the  corresponding  numerical 
solutions.  In  the  tables,  pressure  errors  are  measured  by  Cp  =  l|p  —  Hie  continuous 

norm  of  the  difference.  Because  P  is  piecewise  constant,  first-order  convergence  is  the  best 
that  can  be  expected.  Velocity  errors  are  calculated  separately  for  vertical  and  horizontal 
edges; 


(60) 
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Control  Volume  Mixed 

Finite  Difference 

Grid 

.  .... 

Ilevll 

16x16 

1.54E-4 

9.41E-5 

1.80E-4 

5.65E-4 

6.51E-4 

8.62E-4 

32x32 

8.50E-5 

4.98E-5 

9.85E-5 

3.30E-4 

3.72E-4 

4.98E-4 

64x64 

4.43E-5 

2.51E-5 

5.09E-5 

1.90E-4 

2.09E-4 

2.82E-4 

128x128 

1.90E-5 

1.06E-5 

2.17E-5 

1.08E-4 

1.17E-4 

1.59E-4 

256x256 

— 

— 

— 

6.35E-5 

6.69E-5 

9.22E-5 

Table  1:  Comparison  of  methods  for  uniform  grids  and  variable  permeability  -  entire  domain 


e 


—  V)  '  Uy  ds 
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(61) 


Then  ||ev||  =  (e^  +  is  equivalent  to  a  discrete  H{div)  norm  of  the  vector  velocity 

error  (which  incorporates  the  norms  of  {v  {v-V)y,  and  div(v  -  V),  the  last  of 

which  is  zero  by  the  local  conservation  property  of  the  mixed  method). 

Problem  1.  We  first  compare  the  accuracy  of  the  control-volume  mixed  method  with 
that  of  block-centered  finite  differences  on  rectangular  grids.  Problem  1  subdivides  the 
domain  Q  =  [— 1, 1]^  into  four  quadrants  and  assigns  a  different  value  of  A  (a  scalar)  to  each: 
0.01  for  X  >  0,  y  >  0;  0.05  for  x  <  0,  y  >  0;  10  for  x  <  0,  y  <  0;  33.33  for  a:  >  0,  2/  <  0. 
The  source  term  was  zero  over  the  entire  domain  and  the  boundary  conditions  specified  the 
normal  component  of  the  velocity  on  dQ  as  follows: 


2000/1005  x’rz:-l,y<0, 
10/1005  x  =  -l,y>0, 
6666/3334  x  =  l,y<0, 
2/3334  ar^l,i/>0, 

0  y  =  ±l. 


These  boundary  conditions  specify  both  the  total  flux  in  at  the  left  boundary  and  out  at  the 
right  boundary  to  be  equal  to  2,  and  no  flow  at  the  top  and  bottom  boundaries.  Note  that 
making  the  a^-component  of  the  velocity  on  the  the  left  and  right  boundaries  proportional 
to  the  permeability  avoids  singularities  in  the  velocities  at  (il,0);  however  there  is  still  a 
singularity  at  the  origin.  Uniform  grids  from  16  x  16  to  256  x  256  were  used.  In  the  absence 
of  an  analytical  solution,  the  256  x  256  mixed  solution  was  used  for  comparison.  Velocity 
errors  for  the  two  methods  are  given  in  Table  1.  The  superiority  of  the  mixed  method  is 
evident;  the  accuracy  of  its  velocities  on  a  32  x  32  grid  is  similar  to  finite  differences  on  a 
256  X  256  grid.  This  is  not  an  artifact  of  the  use  of  the  256  x  256  mixed  solution  as  “exact/’ 
since  the  128  x  128  and  256  x  256  finite  difference  solutions  differ  from  each  other  by  an 
order  of  magnitude  more  than  the  corresponding  mixed  solutions  differ  from  each  other. 

This  is  the  one  instance  in  which  the  mixed  velocities  do  not  show  second-order  conver¬ 
gence,  the  reason  being  that  the  true  solution  has  a  singularity  at  the  origin.  Excluding  from 
the  summation  in  Eq.  (60)  those  edges  that  lie  inside  (-|,  i)^  (and  similarly  in 
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Control  Volume  Mixed 

Finite  Difference 

Grid 

Ilevll 

Ilevll 

16x16 

1.28E-4 

9.12E-5 

1.57E-4 

5.57E-4 

4.95E-4 

7.45E-4 

32x32 

3.69E-5 

2.78E-5 

4.62E-5 

2.61E-4 

2.38E-4 

3.53E-4 

64x64 

9.56E-6 

7.46E-6 

1.21E-5 

1.29E-4 

1.19E-4 

1.75E-4 

128x128 

2.06E-6 

1.64E-6 

2.63E-6 

6.42E-5 

5.99E-5 

8.78E-5 

256x256 

— 

— 

— 

3.22E-5 

3.02E-5 

4.42E-5 

Table  2:  Comparison  of  methods  for  uniform  grids  and  variable  permeability  -  away  from 
singularity 


(-1.-1) 

Figure  7:  Macro-blocks  of  distorted  grids 


Eq.  (61)  excluding  those  edges  that  lie  inside  )  we  can  see  the  behavior  of 

the  error  away  from  this  singularity.  Table  2  shows  the  result  when  the  errors  are  calculated 
over  12  \  (— |,  the  domain  minus  a  small  square  centered  at  the  origin.  From  the  table 
we  see  that  the  mixed  method  is  still  superior,  and  that  the  velocity  errors  for  the  mixed 
method  seem  to  indicate  second-order  convergence  away  from  the  singularity  whereas  the 
finite  difference  velocity  errors  appear  to  be  first-order  convergent. 

It  should  be  noted  that  while  this  problem  demonstrates  that  the  mixed  method  can 
produce  more  accurate  velocity  approximations  than  the  finite  difference  method,  the  two 
methods  exhibit  comparable  accuracy  for  the  pressure.  The  remaining  problems  study  the 
accuracy  of  the  control- volume  mixed  method  under  various  conditions. 

Problem  2.  In  this  problem  we  study  the  effect  of  grid  distortion  and  tensor  perme¬ 
ability  on  the  control- volume  mixed  method.  Here  the  distortion  is  based  on  the  angle  (3 
shown  in  Figure  7,  where  finer  grids  are  obtained  by  refining  the  four  macro-blocks  along 
bilinear  coordinate  lines.  The  coefficient  A,  a  constant  anisotropic  tensor,  is  given  by 


COS  6  sin  0 

■  1 

0 

cos  0  —  sin  9 

—  sin  0  cos  0 

0 

0.01 

sin  0  cos  0 

21 


o 

O 

il 

o 

O 

oo 

11 

/?  =  80",6'  =  45° 

Grid 

€p 

INvIl 

€p 

llevll 

16x16 

2.506E-1 

2.789E-3 

2.256E+0 

3.337E-1 

32x32 

1.266E-1 

7.414E-4 

6.332E-1 

1.174E-1 

64x64 

6.353E-2 

1.885E-4 

1.736E-1 

3.355E-2 

128x128 

3.189E-2 

4.734E-5 

5.189E-2 

8.740E-3 

o 

O 

II 

o 

o 

11 

/?  =  60",  =  45° 

Grid 

Cp 

llevll 

Cp 

llevll 

16x16 

2.611E-1 

4.400E-3 

3.081E+0 

5.345E-1 

32x32 

1.334E-1 

1.230E-3 

9.405E-1 

2.338E-1 

64x64 

6.970E-2 

3.225E-4 

2.675E-1 

7.644E-2 

128x128 

3.993E-2 

8.208E-5 

7.817E-2 

2.101E-2 

Table  3:  Mixed  method  accuracy  for  distorted  grids,  constant  tensor  permeability 


where  0  is  the  angle  between  the  coordinate  axes  and  the  principal  directions  of  permeability. 
No-flow  boundary  conditions  were  used  and  the  source  term  was  chosen  to  yield  an  exact 
solution 

p{x,y)  =  cos cos (2 Try).  (63) 

Table  3  presents  results  for  the  extreme  values  of  9  and  several  values  of  /3  as  in  Figure  7. 
The  second-order  convergence  of  the  velocity  is  clear,  though  with  9  =  and  (5  =  60^ 
(serious  anisotropy  and  distortion)  the  asymptotic  regime  is  not  reached  until  the  grid  is 
quite  fine.  We  note  that  the  accuracy  degrades  slightly  as  the  grid  becomes  more  distorted, 
which  is  no  surprise.  The  degradation  of  accuracy  due  to  lack  of  alignment  of  anisotropy 
with  the  grid  is  much  more  severe.  This  effect  is  greater  in  this  test  problem  than  would  be 
expected  in  practice,  since  the  100:1  anisotropy  ratio  exceeds  that  of  typical  porous  media, 
and  a  modeler  would  attempt  to  avoid  the  worst  case  of  ^  =  45^. 

Problem  3.  The  last  problem  considers  a  distorted  grid,  based  on  the  macro-blocks  in 
Figure  8,  along  with  variable  tensor  permeability.  The  tensor  is  given  on  Regions  I,  II,  and 
III  in  Figure  8  by 


A/ 


1/4  1/4  ■ 
1/4  4  J  ’ 


A//  = 


2  1 
1  1  ’ 


A///  = 


2 

1/2 


1/2  ■ 
1/2  J 


(64) 


Boundary  conditions  and  the  source  term  are  specified  so  that  the  exact  solution  is 


Pi{x,y)  =  x'^  -C, 


Pn{x,y)  =  — 


-c 


Piii(x,y)  =  -y^  -  C, 


(65) 


where  C  is  chosen  to  make  the  integral  of  p  vanish.  The  pressure  and  flux  are  continuous 
at  interfaces.  Table  4  reports  the  results,  and  again  we  see  second-order  convergence  in  the 
velocity. 
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5  Conclusions 


The  control-volume  mixed  finite  element  method  provides  a  simple,  systematic,  easily  im¬ 
plemented  means  of  obtaining  accurate,  locally  conservative  velocities  on  irregular  block- 
centered  grids,  allowing  for  effects  of  anisotropy  and  heterogeneity.  Degrees  of  freedom 
consist  of  block  pressures  and  edge  (two  dimensions)  or  face  (three  dimensions)  fluxes,  with 
no  additional  complexities  such  as  Lagrange  multipliers,  so  that  the  method  is  strongly  anal¬ 
ogous  to  block-centered  finite  differences  from  a  modeler’s  point  of  view.  Velocities  obey 
a  rigorously  derived  discrete  Darcy  law  and  exhibit  second-order  convergence  as  long  as 
the  exact  solution  has  no  singularities.  Heterogeneities  and  reasonable  distortions  have  only 
mild  effects  on  the  accuracy  of  the  method.  Severe  anisotropy  that  is  strongly  oblique  to  the 
coordinates  leads  to  significant  increases  in  velocity  errors,  though  they  are  still  second-order 
convergent. 

The  numerical  results  in  Section  4  were  obtained  using  a  multilevel  solver.  The  multilevel 
solver  for  the  discrete  control- volume  mixed  finite  element  equations  is  comparable  in  cost 
to  a  typical  finite  difference  solver.  This  multilevel  solver  is  discussed  in  [16]  and  will  be  the 
subject  of  a  future  paper. 
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